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In the present Letter we use the Direct Numerical Simulation (DNS) of the Navier-Stokes equation 
for a two-phase flow (water and air) to study the dynamics of the modulational instability of free 
surface waves and its contribution to the interaction between ocean and atmosphere. If the steepness 
of the initial wave is large enough, we observe a wave breaking and the formation of large scale dipole 
structures in the air. Because of the multiple steepening and breaking of the waves under unstable 
wave packets, a train of dipoles is released and propagate in the atmosphere at a height comparable 
with the wave length. The amount of energy dissipated by the breaker in water and air is considered 
and, contrary to expectations, we observe that the energy dissipation in air is larger than the one in 
the water. Possible consequences on the wave modelling and on the exchange of aerosols and gases 
between air and water are discussed. 



PACS numbers: Valid PACS appear here 

The modulational instability, also known as the Ben- 
jamin Feir instability, is a well known universal phe- 
nomenon that takes place in many different fields of 
physics such as surface gravity waves, plasma physics, 
nonlinear optics (see the recent historical review pQ). The 
basic idea consists in that a sufficiently steep sinusoidal 
wave may become unstable if perturbed by a long enough 
perturbation. It is a threshold mechanism, therefore, for 
example, for surface gravity waves in infinite water depth 
a wave is unstable if 2y/2koAo > Afc/fco, where fen is the 
wave number of the sinusoidal wave (carrier wave), Ak 
is the wave number of the perturbation and A is the 
amplitude of the initial wave. 

The modulational instability has been discovered in 
the sixties and recently it has received again attention be- 
cause it has been recognized as a possible mechanism of 
formation of the rogue waves [3] . The standard math- 
ematical tool used to describe such physical phenomena 
is the Nonlinear Schrordinger equation (NLS), which is 
a weakly nonlinear, narrow band approximation of some 
primitive equation of motion. The beauty of such equa- 
tion is that it is integrable and many analytical solutions 
can be written explicitly. For example, breather solutions 
[UG3 have been considered as prototypes of rogue waves; 
they have been observed in controlled experiments both 
in surface gravity waves and in nonlinear optics [BHH]- 

Concerning the ocean waves, the studies on the mod- 
ulational instability have been concentrated on the NLS 
dynamics and only more recently numerical computation 
of non-viscous, potential, fully nonlinear equations have 
been considered jTDl E]. However, so far none of the 
aforementioned literature has ever considered the effect 
of the modulation instability on the fluid above the free 
surface. As far as we know this is the first attempt in 
which the dynamics of air on water during the modula- 



tion process is investigated. This has been possible by 
simulating the Navier-Stokes equation for a two phase 
flow. This approach allows us to investigate conditions 
which are beyond the formal applicability of the NLS 
equation: for example, it is well known that if the initial 
wave steepness is large enough, the NLS equation is not 
able to describe the dynamics because the breaking of the 
wave takes place jTSJ, [T3] . For steep waves and particu- 
larly those close to the breaking onset, vorticity is gen- 
erated by viscous effects and by the topological change 
of the interface in case of bubble entrainment processes 
(this dynamics cannot be described by potential models) . 
Breaking of surface waves, as an oceanic phenomenon 
[14], is important across a very broad range of applica- 
tions related to wave dynamics, atmospheric boundary 
layer, air-sea-interactions, upper ocean turbulence mix- 
ing, with respective connections to the large-scale pro- 
cesses including ocean circulation, weather and climate 
15 . Modulational instability and breaking has become 
also relevant in engineering applications [16] such as naval 
architecture, structural design of offshore developments, 
marine transportation, navigation, among many others. 

In the present work, the two-fluid flow of air and 
water is approximated as that of a single incompress- 
ible fluid with density and viscosity smoothly varying 
across the interface. The continuity and momentum 
equations (Navier-Stokes) are solved in generalized co- 
ordinates [T71 HH] . The variation of the fluid properties 
and the surface tension forces is spread across a small 
neighborhood of the interface. The interface between air 
and water is captured as the zero level-set of a signed 
distance from the interface d(x,t) which, at t = 0, is 
initialized by assuming d > in water, d < in air. 
Physical fluid properties are assumed to be related to d 
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by the equation: 

f(d) = f a + (f W -fa)H S (d) (1) 

where Hg(d) is a smooth step function and the parameter 
6 is chosen so that the density and viscosity jumps are 
spread across some grid cells |19) . The distance function 
is advected in time with the flow as a non-diffusive scalar 
by using the equation 

dd 

-r+u-Vd = 0, (2) 
at 

and the interface is located as the d = level. Hence, 
the distance is reinitialized. 

In our simulations we consider the standard modula- 
tional instability process as the one produced for example 
in the experimental work in |20j . The initial surface el- 
evation is characterized by a perturbed sinusoidal free 
surface elevation of the form: 

ri(x, t = 0) = A Q cos(k x) + Ai(cos(k + x) + cos(k~x)), 

(3) 

where k is the wave number of the carrier wave, k — 
ko ± Ak with Afc the wavenumber of the perturbation. 
The simulations presented hereafter are characterized by 
a eo = koAo that is varied from 0.1 to 0.18, with a step 
0.02. The sideband components are placed at Ak = k /5 
and their amplitude is A\ — 0.1-Ao. It is worth notic- 
ing that the conditions are essentially similar to those 
used in [11] [2Tj and corresponds to the early stages of an 
Akhmediev breather [2"2"] . 

Because of the typical time scale of the modulational 
instability is of the order of 100 periods, the Navier- 
Stokes simulation is expensive and for the initial devel- 
opment of the instability a standard potential code is 
used. The Navier-Stokes simulation is then initialized 
with solution from the potential flow. For convenience, 
results are presented in dimensional form. Simulations 
are carried out for carrier wave of wavelength Ao = 0.60 
m, with g — 9.81 m s . The computational domain 
which spans horizontally from x — -1.5 m to x = 1.5 
m and vertically from y= -2 m up to 0.6 m above the 
still water level. The domain is uniformly discretized in 
the horizontal direction with Aa; = 1/1024 m. Verti- 
cally, the grid spacing is uniform, and equal to Ax, from 
y = -0.15 to 0.15 m whereas it grows geometrically by 
a factor a— 1.03 towards the upper and lower bound- 
aries. This gives a total of 3072 x 672 grid cells. The 
total thickness of the transition region is 0.01 m, so that 
the density jump is spread across about 10 grid cells. 
Although neglected in the potential flow model, surface 
tension effects are considered in the Navier-Stokes solu- 
tion for which the surface tension coefficient is assumed 
to be that in standard conditions a = 0.073 N m . In 
the two-fluid modelling, the densities of air and water are 
the standard ones, g w = 1000 kg m -3 and g a = 1-25 kg 
m~ 3 . The values of the dynamic viscosities in water and 
air are /j, w = 10 -3 kg m -1 s _1 and \i a = 1.8 10~ 5 kg m _1 
s _1 , respectively. 
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FIG. 1: Breaking event for initial steepness eo= 0.18. The 
color corresponds to the vorticity field: green corresponds to 
zero vorticity, blue to -20 sec -1 and red to 20 sec -1 . 



The process observed in the simulation corresponds 
to the standard modulational instability (exponential 
growth of the side bands) up to the point where the wave 
group reaches its strongly nonlinear regime and eventu- 
ally wave breaking is observed. We first concentrate our 
attention to the breaking event: in figure [l] an example 
for steepness eo = koAo— 0.18 of the simulation of the 
wave breaking resulting from the modulational instabil- 
ity is shown; the formation of the first jet is observed 
which entraps the air. The jet then bounces on the free 
surface creating a second air bubble. Some droplets of 
water in air are also visible; a small amount of vorticity 
is also released beneath the surface. 

It is of particular interest to focus the attention on 
the dynamics of the air flow. In figure [2] we show a se- 
quence of snap-shots of the water and air domain where 
the formation and detachment of a dipole structure is ob- 
served (the simulation is performed for initial steepness 
of e = 0.16): the rather fast steepening of the wave pro- 
file and the wave breaking causes the air flow to separate 
from the crest giving a rise to large, positive, vorticity 
structure (positive vorticity is in red). The interaction 
of this vortex structure with the free surface leads to 
the formation of a secondary vorticity structure of oppo- 
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FIG. 2: Sequences of formation of a dipole in air as a con- 
sequence of the wave breaking. Results refer to an initial 
steepness eg — 0.16 
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FIG. 3: Vorticity field in a portion of the computational do- 
main for initial steepness eo = 0.18. The scale of vorticity 
is as described in the label of figure [l] The supplemental ma- 
terial contains the animation of the simulation (file: anima- 
tionl8.avi). 
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FIG. 4: Dissipation nondimensionalized with period T and 
total energy in the water, Eg a t = as a function of nondi- 
mensional time for water (solid line) and air (dashed line). 
The initial steepness is 0.12 



site sign, which eventually detaches from the free surface 
and forms a dipole which moves under the self-induced 
velocity (such phenomenon is not observed if the wave 
breaking process does not take place: we have performed 
a numerical simulation with eo = 0.1 which does not lead 
to a breaking event and no evident dipole structures have 
been observed). 

Because the group velocity is half the phase velocity, 
each single wave that passes below the group (at its maxi- 
mum height), breaks. The result is that a series of dipoles 
are released into the atmosphere as shown in figure [3] 
Vortices of various sizes and dipoles are clearly observ- 
able in the domain. Two things should be noted: i) the 
height of the highest dipoles is of the order of the wave- 
length; ii) large amount of vorticity is observed in the 
air and not in the water. One major question to be 
answered, especially in the spirit of modelling the dissi- 
pation term in the wave forecasting models [UJ [21] is the 
amount of energy dissipated during a wave breaking or 
a sequence of breaking events. Therefore, a quantitative 
estimate of the dissipated energy both in air and water 
can be obtained by integrating the viscous stresses over 
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FIG. 5: Dissipation nondimensionalized with period T and 
total energy in the water, Eg a t = as a function of nondi- 
mensional time for water (solid line) and air (dashed line). 
The initial steepness is 0.18 
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FIG. 6: Integrated dissipation for different values of the steep- 
ness in water (w) and air (a). 

the air and water domain: 

e% s3 (t)= 2^— Idardy, (4) 
Jd>8 ax j 

£disa(t)= 2/zey— -i dx dy, (5) 

Jd<-8 ax j 

where is the symmetric part of the strain tensor. In 
figure [4] and [5] we show the dissipation function normal- 
ized by the initial energy of the water and wave period, 
T, as a function of time, nodimensionalized by T, for 
a simulation with steepness eo = 0.12 and eo = 0.18, re- 
spectivelly; the origin of the time axis is set to the time at 
which the Navier Stokes simulations takes over the simu- 
lation with the potential code. The occurrence of spikes 
in the energy content in the air indicates that an en- 
ergy fraction is transferred to the air and is successively 
dissipated by the viscous stresses. The multiple peaks 
correspond to the the multiple breaking that occurs dur- 
ing the modulational instability process (as mentioned 
before, the group is slower than the phase velocity and 
each wave breaks as it passes below it). Similar plots are 
observed for larger steepness. 

In order to show how relevant is the energy dissipa- 
tion in the air in comparison to the corresponding one in 
water, we consider the following integrated quantity: 

Ediss(t) = / e diss (t')dt' '. (6) 
Jo 

The integral is considered for both for water and air. His- 
tories of the integrals in time of the viscous dissipation 



terms in the two media are shown in figure [6j It is rather 
interesting to see that, in nondimensional form, the so- 
lutions for the four different cases almost overlap, with 
a total dissipation in the air about three times that in 
water. 

Discussion - Modulated waves of different initial steep- 
nesses have been analyzed using the two phase flow NS 
equation. For initial steepness larger or equal to 0.12, 
multiple breaking events have been observed. Spray is 
the natural consequence of the wave breaking. Droplets 
of water are thrown in the air; some of these particles 
are so small (aerosols) that they can remain in the air for 
a very long time, forming condensation nuclei for clouds 
and affecting incoming solar radiation and therefore rele- 
vant for climate modelling. Vortices, observed in our sim- 
ulations, can in principle transport aerosols (not resolved 
in our simulations) up to the height of the wave lengths 
(this can be even underestimated because of the presence 
of the solid boundary at the top of the computational do- 
main). Quantitative analyses of the energy contents and 
of the viscous dissipation term have been provided. Even 
bearing in mind the limits of the numerical scheme, the 
results indicate that, due to the highly rotational flow in 
air, the dissipation of the energy is mostly concentrated 
in the air side. We stress that it is a common practice 
to estimate the energy loss due to a wave breaking by 
looking at the amount of energy dissipate in the water, 
see for example [25]. Such measurement are the bases 
for the construction of the dissipation function in the 
operational wave forecasting modelling [53]. If the mod- 
ulational instability is responsible of the wave breaking 
in the ocean, then such results would clearly lead to an 
underestimation of the the total energy dissipated. The 
present work represents the first attempt to study the 
modulational instability starting from the Navier-Stokes 
equation and its effects on water and air. Clearly, it has 
a number of limitation due to the heaviness of the com- 
putation. Probably the most important one is that we 
have assumed that waves are long crested and the fluid 
domain is 2-dimensional. In reality, we expect that 3D 
effects can take place and dipoles and vortices can be- 
come unstable. We also underline that our simulations 
correspond to the propagation of waves without the pres- 
ence of external wind: what would be the consequences 
of a turbulent wind on the generation of vorticity during 
breaking event is under investigation. 

Acknowledgments M.O. has been funded by EU, 
project EXTREME SEAS (SCP8-GA-2009-234175) and 
by ONR grant N000141010991. Dr. Proment and Dr. 
Giulinico are acknowledged for discussion. The work 
by A.I. has been done within the RITMARE research 
project. 



[1] V. Zakharov and L. Ostrovsky, Physica D: Nonlinear [2] M. Onorato, A. Osborne, M. Serio, and S. Bertone, Phys. 
Phenomena 238, 540 (2009). Rev. Lett. 86, 5831 (2001). 



5 



[3] P. A. E. M. Janssen, J. Phys. Ocean. 33, 863 (2003). 
[4] A. Osborne, M. Onorato, and M. Serio, Phys. Lett. A 

275, 386 (2000). 
[5] K. B. Dysthe and K. Trulsen, Physica Scripta T82, 48 

(1999). 

[6] A. Chabchoub, N. Hoffmann, and N. Akhmediev, Phys- 
ical Review Letters 106, 204502 (2011). 

[7] G. Clauss, M. Klein, and M. Onorato (Proceedings of 
the ASME 2011 30th International Conference on Ocean, 
Offshore and Arctic Engineering OMAE2011, June 1924, 
2011, Rotterdam, The Netherlands., 2011). 

[8] B. Kibler, J. Fatome, C. Finot, G. Millot, F. Dias, 
G. Genty, N. Akhmediev, and J. Dudley, Nature Physics 
6, 790 (2010), ISSN 1745-2473. 

[9] B. Kibler, J. Fatome, C. Finot, G. Millot, G. Genty, 
B. Wetzel, N. Akhmediev, F. Dias, and J. Dudley, Scien- 
tific reports 2 (2012). 
[10] A. Dyachenko and V. Zakharov, JETP letters 88, 307 
(2008). 

[11] A. Babanin, D. Chalikov, I. Young, and I. Savelyev, Geo- 

phys. Res. Lett 34, L07605 (2007). 
[12] K. Henderson, D. Peregrine, and J. Dold, Wave Motion 

29, 341 (1999). 

[13] A. Babanin, T. Waseda, T. Kinoshita, and A. Toffoli, 
Journal of Physical Oceanography 41, 145 (2011). 

[14] A. Babanin, Breaking and dissipation of ocean surface 
waves (Cambridge Univ Pr, 2011). 

[15] L. Cavaleri, B. Fox-Kemper, and M. Hemer, Bulletin of 
the American Meteorological Society, Accepted for pub- 



lication (2012). 

[16] G. F. Clauss, M. Klein, M. Dudek, and M. Ono- 
rato (Proceedings of the ASME 2012 31th International 
Conference on Ocean, Offshore and Arctic Engineering 
OMAE2012, July 16, 2012, Rio de Janeiro, Brazil, 2012). 

[17] A. Iafrati, Journal of Fluid Mechanics 622, 371 (2009). 

[18] A. Iafrati, Journal of Geophysical Research 116, C07024 
(2011). 

[19] A. Iafrati and E. Campana, Journal of Fluid Mechanics 

529, 311 (2005). 
[20] M. P. Tulin and T. Waseda, J. Fluid Mech. 378, 197 

(1999). 

[21] J. Song and M. Banner, J. Phys. Ocean. 32, 2541 (2002). 
[22] N. Akhmediev, V. Eleonskii, and N. Kulagin, Theoretical 

and Mathematical Physics 72, 809 (1987), ISSN 0040- 

5779. 

[23] G. Komen, L. Cavaleri, M. Donelan, K. Hasselmann, 
H. Hasselmann, and P. Janssen, Dynamics and modeling 
of ocean waves (Cambridge University Press, Cambridge, 
1994). 

[24] L. Cavaleri, J. Alves, F. Ardhuin, A. Babanin, M. Ban- 
ner, K. Belibassakis, M. Benoit, M. Donelan, J. Groe- 
neweg, T. Herbers, et al., Progress in Oceanography 75, 
603 (2007). 

[25] E. Terray, M. Donelan, Y. Agrawal, W. Drennan, 
K. Kahma, A. Williams, P. Hwang, and S. Kitaigorod- 
skii, Journal of Physical Oceanography 26, 792 (1996). 



